/* This function accepts the "nxpt" number of grid-point co-ordinates "x" to determine the
   co-efficine matrix. It fills the structure "cm" with the co-efficients corresponsing to
   "ndiscr" order central-difference discretization of "nderv"-th derivative.

   * If xderv = 0, the derivative is evaluated at the points given by "x".

   * If xderv != 0 the derivative is evaluated at the points given by "xderv".
      -->> Size of the "xderv" array has to be "nxpt-1".
      -->> If xderv[i] = (x[i] + x[i+1])/2, the most accurate derivativs are calculated.
           As it deviates from the above condition, the derivative becomes more and more errorineous.
           This is a better idea to adjust the array "x" so that "xderv" satisfy the above condition.

   * Size of the "tempmtrx" is [nxpt][nxpt]. This is temporary space to speed-up the evaluation.
     If it is not provided, it is allocated dynamically and distroyed after the use.

   * The memory allocated to the elements of matrix "cm" needs to be distroyed after use in the calling function.

   * The p-th row of the cm->mtrx corresponds to istart[p]~iend[p]-th elements of the u[0~nxpt] vector.
     Here, p = 0, 1, 2, ...., cm->nrows.

   * Range of the p-th row of the cm->mtrx: 0 ~ iend[p]-istart[p].

   * Size of the p-th row of the cm->mtrx: iend[p]-istart[p]+1.

   * The derivative can be calculated as: derv[q] = cm->mtrx[q][0 ~ iend[q]-istart[q]]*u[istart[q]~iend[q]].
     Here, q = 0, 1, 2, ..., cm->nrows. */

# include "scene.h"
# include <stdlib.h>

extern short cmtrx_CD(double *x, unsigned short nxpt, double *xderv, unsigned short nxdervpt,\
 unsigned short ndiscr, double **mtrx);

extern short cmtrx_make(double **mtrx, unsigned short nrows, unsigned short ncols, struct cmtrx *cm);

extern short cmtrx_nthderv_CD(double *x, unsigned short nxpt, double *xderv, struct cmtrx *cm, double **tempmtrx[3])
  {
   unsigned short i, j, k, iderv, nxpt2, nxdervpt2, ndiscr;
   double *x2, *xderv2, ***_cmtrx, *xmid;
   short ecode, nderv;

   nderv = cm->nderv;
   ndiscr = cm->napprox;

/* The order of discretization has to be even. */
   if(ndiscr%2 != 0) return -10;

/* Memory allocation for temporary use. */
   if(tempmtrx != (double ***)0) _cmtrx = tempmtrx;
   else{
        _cmtrx = (double ***)malloc(sizeof(double **)*3);
        for(i = 0; i < 3; i++)
           {
            _cmtrx[i] = (double **)malloc(sizeof(double *)*nxpt);
            for(j = 0; j < nxpt; j++) _cmtrx[i][j] = (double *)malloc(sizeof(double)*nxpt);
           }
       }

   xmid = (double *)0;
   if(xderv != (double *)0) xmid = xderv;
   else if(nderv > 1)
       {
        xmid = (double *)malloc(sizeof(double)*(nxpt-1));
        for(i = 0; i <= nxpt-2; i++) xmid[i] = 0.5*(x[i]+x[i+1]);
       }

/* Start with an unit matrix. */
   for(i = 0; i < nxpt; i++) for(j = 0; j < nxpt; j++) _cmtrx[1][i][j] = (double)((i == j) ? 1:0);

   for(iderv = 1; iderv <= nderv; iderv++)
     {
   /* This determines the location at which iderv-th derivative is to be avaluated. */
      if( (iderv + nderv + (xderv == (double *)0 ? 0:1))%2 == 1)
          {
        /* The iderv-th derivative will be evaluated at points given
           by "xmid" after using values at "x". */
           x2 = x;
           nxpt2 = nxpt;
           xderv2 = xmid;
           nxdervpt2 = nxpt-1;
          }
      else{
        /* The iderv-th derivative will be evaluated at points given by "x"
           after using values at "xmid". */
           x2 = xmid;
           nxpt2 = nxpt-1;
           xderv2 = x;
           nxdervpt2 = nxpt;
           if(iderv == 1)
              {
               x2 = x;
               nxpt2 = nxpt;
              }
          }

   /* The co-efficient matrix of the 1st-derivative is returned the 2D-array "_cmtrx[0]". */
      if( (ecode = cmtrx_CD(x2, nxpt2, xderv2, nxdervpt2, ndiscr, _cmtrx[0])) < 0) return ecode;

   /* Multiplying the newly-obtained matrix with the last coefficient-matrix for evaluation of
      more than 1st-order derivatives. It uses the following idea:
      {u'}  = [A1]{u}
      {u''} = [A2]{u'} = [A2][A1][{u} = [B1]{u}. etc.

      This is the final required matrix. */
      for(i = 0; i < nxdervpt2; i++) for(j = 0; j < nxpt; j++)
      for(_cmtrx[2][i][j] = 0, k = 0; k < nxpt2; k++) _cmtrx[2][i][j] += _cmtrx[0][i][k]*_cmtrx[1][k][j];

      for(i = 0; i < nxdervpt2; i++) for(j = 0; j < nxpt; j++) _cmtrx[1][i][j] = _cmtrx[2][i][j];
     }

   if(xderv == (double *)0) if(xmid != (double *)0) free(xmid);

/* This function filters out all the zero-elements and keeps only non-zero elements to
   save the memory and to reduce the processing time during multiplication etc. */
   cmtrx_make(_cmtrx[2], nxdervpt2, nxpt, cm);

/* Free the temporary allocated memory. */
   if(tempmtrx == (double ***)0)
     {
      for(i = 0; i < 3; i++)
        {
         for(j = 0; j < nxpt; j++) free(_cmtrx[i][j]);
         free(_cmtrx[i]);
        }
      free(_cmtrx);
     }

   return 1;
  }
